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Abstract 

The photon scattering properties of aperiodic nano-scale dielectric structures can be tailored 
to closely match a desired response by using adaptive algorithms for device design. We show 
that broken symmetry of aperiodic designs provides access to device functions not available to 
conventional periodic photonic crystal structures. 
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A. Introduction 



The spatial arrangement of nano-scale dielectric scattering centers embedded in an oth- 
erwise uniform medium can strongly influence propagation of an incident electromagnetic 
wave. Exploiting this fact, we iteratively solve an inverse problem Q to find a spatial ar- 
rangement of identical non-overlapping scattering centers that closely matches a desired, or 
target, response. Of course, the efficiency of adaptive algorithms used to find solutions may 
become an increasingly critical issue as the number of scattering centers, N, increases. How- 
ever, even for modest values of N, this method holds the promise of creating nano-photonic 
device designs that outperform conventional approaches based on spatially periodic photonic 
crystal (PC) structures. 

The configurations considered by us consist of either lossless dielectric rods in air or 
circular holes in a dielectric similar to the majority of quasi two dimensional (2D) PCs 
reported in the literature by, for example, [2, 3J. To confirm the validity of our 2D simulations 
we compare the optimized configuration with full three dimensional (3D) simulations. 

B. Calculation of the scattered field 

The optimization process is an iterative procedure, in which the scattered field from a 
trial configuration of cylindrical rods or circular holes is compared with that of the target. 
We note that computation of the scattered field is sometimes referred to as the forward 
problem. 

The electromagnetic field solver used is based on the analytical solution of the Helmholtz 
equation by separation of variables in polar coordinates. A typical problem we consider is 
a set of N long, parallel, lossless circular dielectric rods distributed in a uniform medium 
and illuminated by an electromagnetic wave perpendicular to the axis of the cylinders. The 
natural geometry of the system led us to use a 2D electromagnetic field solver. When 
analyzing scattering from one cylinder, the solution, expressed as a Fourier-Bessel series, is 
found immediately by imposing the continuity of the electric and magnetic field components 
at the rod surface. However, when studying scattering from two or more cylinders, multiple 
scattering gives an additional linear system that has to be solved in order to find the Fourier- 
Bessel coefficients j^. For a given number of Bessel functions used this linear system has a 
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reduced form which is conveniently described using the scattering matrix method 

0]. The input wave can, in principle, be of arbitrary shape so long as it is expressed as 
a Fourier-Bessel series. In our simulations the input beam is a gaussian and both TE and 
TM polarizations are considered. The calculated Fourier-Bessel coefficients for the gaussian 



beam conform 
in Appendix lAl 
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111 ]. Additional details on the electromagnetic solver used may be found 



C. Optimization 

For the results presented in this paper, the optimization method used is the guided 
random walk. The positions of individual cylinders are randomly changed by a small amount, 
the modified distribution is simulated, and if the result is closer to the target function, such 
that the error defined in Appendix B is decreased, then the new configuration is accepted, 
otherwise it is rejected. The actual implementation of the adaptive algorithm also includes 
other types of collective motion such as moving more than one cylinder per iteration, and 
moving or rotating all the cylinders. 

In general the target function may be any electromagnetic distribution. A typical target 
function might involve redirecting and reshaping the input beam. More complex situations 
involving optically active materials are not considered in this paper. Because the target 
function can specify the field distribution of a guided mode, it is possible to create mode 
converters to interface between a fiber or ridge waveguide and a PC waveguide. In this paper 
we report on redirecting and reshaping an electromagnetic beam using optimized aperiodic 
nano-structures. 



D. Results 

To illustrate our approach, we considered an input beam of gaussian profile scattered by 
an angle of 45°. The target functions we consider are a top hat distribution of the optical 
intensity with respect to the scattering angle and a cosine squared (cos 2 ) distribution of the 
intensity. We note that a modal field distribution target function requires amplitude and 
phase to be specified. 

As an initial demonstration we chose the top-hat intensity function because it is difficult 
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to achieve (even approximately) in conventional optical systems and serves to show the 
capabilities of our adaptive algorithm. Application of such an intensity distribution includes 
guaranteeing uniform illumination of the active area of a photodetector. On the other hand 
the cos 2 intensity distribution target approximates the transverse spatial mode intensity 
typically found in a waveguide and so could be considered a first step toward design of a 
waveguide coupler. 

1. Top-hat target function 

Optimization of a top-hat intensity distribution target function is performed starting 
from a configuration of N = 56 dielectric rods (represented in Fig. CJa) by the small 
circles), each having refractive index n = 1.5, and diameter d = 0.4 /im. The medium 
surrounding the rods is air and the structure is illuminated by a TM polarized (electric 
field along the z direction) gaussian beam of width 2a = 4 ^m, wavelength A = 1 /im, 
and propagating along the positive x direction (from left to right in Fig. IHa)). The initial 
configuration of the rods and intensity distribution are illustrated in Fig. HJa), where the 
arrows represent (in arbitrary units) the real part of the Poyting vectors. The target function 
window (represented in Fig. ^a) by a missing arc in the 7 /im radius observation circle) 
extends from 30° to 60°. Fig. Q (b) shows normalized intensity (real part of the normal 
component of total field Poynting vectors directed outwards) as a function of angle on a 
radius of 7 fim from the center of the symmetric array. Clearly, for the initial configuration, 
the overlap with the top-hat target function (broken line) is poor. 

Fig. |2f a) shows the spatial distribution of the N = 56 rods with the real part of Poynting 
vectors after 9700 iterations of the adaptive search algorithm. Fig. I2fb) shows the corre- 
sponding angular distribution of the intensity. In Fig. |2(c) the distribution of the electric 
field (relative magnitude, with 1 corresponding to the maximum magnitude in the incident 
gaussian beam) is shown and in Fig. I2fd) the relative error versus number of iterations (the 
errors are normalized with respect to the initial value) are represented. Note that for this 
system with 2 x N positional degrees of freedom, the error is not saturated even after 9,700 
iterations (Fig. Efd)) and the structure can be further optimized by performing additional 
iterations. The error is computed using the method described in Appendix El with exponent 
7 = 2 on the observation circle of radius 7 //m. 



4 



For comparison, in Fig. |3] we show the results of using a PC with the same number of 
dielectric rods, N = 56. Clearly, the number of rods is not sufficient to redirect all the beam 
in the 45° (±15°) direction and the error with respect to the target function is unacceptably 
large. The spatial symmetry of the PC excludes access to the top-hat target function. It 
is only by breaking the symmetry that one may acquire the target. In general, broken 
symmetry enables functionality. 



2. Cosine squared target function 

For a cos 2 target function we chose a different situation in which there are N = 26 lower 
index cylinders (Si02, n = 1.45) embedded in a higher index material (Si, n = 3.5). The 3D 
equivalent of the modelled situation is a Si slab with cylindrical perforations embedded in 
Si02- The incident wave is a TE polarized (magnetic field along the z direction) gaussian 
beam of width 2a = 1.5 /iin and wavelength A = 1.5 /mi. In a similar manner to the 
previous example, the initial configuration and intensity distribution are illustrated in Fig. 
0] The back scattered intensity in the input region is calculated as the difference between 
the incident field Poynting vector and total field Poynting vector. In Fig. we show the 
optimized configuration, Poynting vectors, relative error versus number of iterations, and 
relative magnitude of the magnetic field (with 1 corresponding to the maximum magnitude 
in the original gaussian beam). This time the position of the N = 26 cylinders could be 
changed and so could their diameters d, giving a total number of degrees of freedom that 
is three times the number of cylinders N. The values of the diameters are constrained to 
the range 0.2 < d < 0.5 /mi. The error is computed using the metric from Appendix IB1 with 
exponent 7 = 1 on a 6 /tm radius circle. Notice that due to the small number of cylinders 
the actual number of degrees of freedom is smaller (78) relative to the previous calculations 
with a top-hat target function (112) and the optimization saturates after a comparatively 
small number of iterations. 



3. Computing resources 

Our previous optimization work on ID problems required relatively insignificant 
computational resources. Optimization algorithms using 2D electromagnetic solvers are 
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always more compute intensive than those for ID structures. Interesting simulations in 
3D are even more compute intensive and often require parallel computing. Each of the 
simulations (top-hat and cos 2 target) discussed in this paper took 4 days to complete (9,700 
iterations in the first case and 7,600 of the second) using a Pentium IV processor with a 3 
GHz clock frequency, 533 MHz memory bus, and 1 GB RDRAM. 

The compute time for the forward problem solver is dominated by the solution of a linear 
system with a full matrix of complex numbers and thus is strongly dependent on the number 
of cylinders N considered in the problem and the number of Bessel functions Nj, used. The 
size of the system is N m = N(2Nb + 1) and the solver routine (iteratively corrected LU 
decomposition Q|) time is proportional to the cube of the matrix size (O(iV^J). When 
using an incident plane wave or a cylindrical wave the number of Bessel functions needed is 
very small, however, an appropriately accurate approximation of the gaussian beam shape 
requires a large number of Bessel functions and a corresponding increase in compute time. 



4- Comparison with 3D Simulation 

For practical purposes, a comparison with 3D electromagnetic simulations serves to con- 
firm the accuracy of solutions obtained with the 2D simulator. Also, since 3D simulations 
are much more time consuming than 2D ones, one might adopt the 2D optimization result 
as a starting point for a limited number of 3D iterations which are then used to further 
refine optimization. 

In this paper we compare results of our 2D electromagnetic simulations with 3D Finite 



Integration Technique ( 
CST Microwave Studio 



simulations obtained using a commercially available package, 



In a realistic structure one might anticipate the infinitely long cylindrical hole structure 
used in the cos 2 target simulation to be replaced with Si02 filled holes in a Si slab itself 
embedded in Si02- The input beam might be launched into this slab from a ridge waveguide. 
In our simulations, the wavelength of the light is A = 1.5 /im and the polarization is TE. The 
Si slab is 0.6 fxm thick and the effective index of the fundamental mode of this slab waveguide 
is the same as the index of the material surrounding the cylinders in the 2D simulation. The 
diameter of the holes is d = 0.4 /mi and the mode size of the ridge waveguide (Fig. [HJ) is 
approximately 1.5 fim and so the same as the width of the gaussian beam used for the 2D 
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simulations. The differences between the two simulations (Fig. EJb) for the 3D and Fig. 
EJd) for the 2D simulation results) are primarily due to the nonuniformity of the field and 
structure in the z direction but also because of the difference between the gaussian beam 
and the ridge mode and the fundamently 3D discontinuity between the ridge waveguide and 
the slab. These differences translate into a 30% decrease in peak intensity and 25% decrease 
in the total power directed in the desired direction for the 3D simulation compared to the 
2D one (Fig. El (c)). 

5. Sensitivity analysis 

We analyzed the influence of small changes in the wavelength of the incident beam on 
error. We took the optimized configuration from the cos 2 target example (Fig. Efa)) and 
simulated the effect of changing the frequency of the input light. The relative deviation from 
the minimum error is shown in Fig. as a function of the electromagnetic wave frequency 
shift. Frequency variations of Af = 200 GHz only change the error by 0.3% of its minimum 
value. Thus an optical beam centered at wavelength A = 1.5 /im modulated at very high 
speed will behave essentially as the simulated monochromatic wave at A = 1.5 fim (200 
THz). Even a 1 THz deviation in frequency changes the error by only 6 - 7 %. 

Another important issue is the sensitivity of the aperiodic nano-photonic design to vari- 
ations and errors introduced by the fabrication process. As an initial study we explored 
the sensitivity of the design to changes in the position of the cylinders. Other parameters 
that are influenced by fabrication processes such as the diameter and detailed shape of the 
cylinders, index of refraction, direction of the input ridge waveguide with respect to the 
cylinders, etc., are not analyzed in this paper. 

To estimate the sensitivity to position, for the same cos 2 target function example, we use 
the following method. First, one by one the cylinders are displaced by a small fixed distance 
A in the positive and negative x and y directions and the change in error is evaluated for 
each cylinder (the maximum error created out of the 4 displacements). The displacement 
values used are small compared to the initial diameters of the cylinders (d = 400 nm). Five 
values for the displacement A are considered: 10, 20, 30, 40, and 50 nm. Next the 10 
displaced cylinders with the greatest influence on the error are selected. These cylinders, as 
expected, are located in regions associated with high field intensity. 
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The selected cylinders are all individually moved randomly by the same step size in the 
positive and negative x and y directions (each movement has the same 1/4 probability). The 
movement of each cylinder is independent of the movements of the other selected cylinders. 
This way we generate a number of slightly modified configurations. The errors for 10 of 
these modified configurations are computed for each displacement value and a combined 
plot generated (Fig. |HJ). 

We chose this very simplified method for estimating the sensitivity to position because 
a more complete approach would involve independent displacements of each cylinder in 
completely random directions and with variable distances and hence give rise to significant 
computational effort. The plot in Fig. |S1 suggests that 10 nm precision in fabrication 
might be needed to ensure less than 10% increase in minimum error for devices operating 
at wavelength A = 1.5 /im. 



6. Comparison with photonic crystal inspired devices 

In recent years, spatially periodic dielectric structures have been studied and applied 
to both optics and microwaves Q, Q] . It has been shown that introduction of point 
and line defects in PCs can be used to filter, demultiplex, and guide electromagnetic waves 



Il9| . However, there are numerous unresolved design issues with PC inspired devices 
that limit prospects of adoption as practical components. For example, when coupling 
between standard fibers or waveguides and PC waveguides, the back reflection is either 
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221. As 



unacceptably large |20j or requires use of relatively large coupling regions 
another example, it is well known that the necessarily finite size of the PC can have a 
dramatic and detrimental impact on device performance j^J . Solutions to these and similar 
problems are stymied by the limited number of degrees of freedom inherent to PC design. 
Rather than struggle for solutions within the constraints of spatial symmetry imposed by 
PC structures, our approach is based on breaking the underlying spatial symmetries and 
thereby exposing larger numbers of degrees of freedom with which to design and optimize 
nano-photonic devices. 

Our experience so far indicates that optimization of such systems is best achieved using 
numerical adaptive design techniques. This is because solutions, such as that illustrated 
in Fig. 03 have such a high degree of broken symmetry it is unlikely analytic methods or 
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conventional intuition could usefully be applied. 

7. Conclusions 

In this paper we have shown that aperiodic nano-photonic dielectric structures designed 
using adaptive algorithms can be tailored to closely match desired electromagnetic trans- 
mission and scattering properties. It is the broken symmetry of the structure that allows 
more degrees of freedom and the possibility of better optimization compared to symmetric 
photonic crystal structures. 

In general the frequency response and spatial response of this system can have very 
complicated forms. It is the large number of degrees of freedom that allow us to tailor the 
response to the desired target. The configuration space is even more complex if materials 
with optical loss or gain are considered. 
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APPENDIX A: ELECTROMAGNETIC SOLVER 

For brevity we will only discuss TM polarized electromagnetic waves. For the TE case 
the equations are similar with H replacing E. The total field is written as the sum of the 
incident field E inc and the field E sc scattered from the cylinders E = E inc + J^iLi^L.- The 
actual incident field on a cylinder labled with index j is E{ nc = E inc + !» m- We 

are interested in solving the Helmholtz equation for the total field, V 2 E + k 2 E = 0, where 
k = ko in the region outside the cylinders and k — k\ inside the cylinders. The method 
used to solve this equation is separation of variables in polar coordinates. Of course, the 
solution of the homogeneous Helmholtz equation in polar coordinates is well known. All 
field quantities may be written in the form of Fourier-Bessel series with the coefficients a m 
and p m determined from the boundary conditions. Hence, E = J2m=-oo a m^m{kp)e irnd + 
Ylm=-oo Pm,Zm(kp)e irn6 , where Z m and Z m are two conjugate cylindrical functions. These 
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functions are either the first order Bessel functions J m and Y m or the second order Bessel 
functions (Hankel functions) Hm and Hm ■ The pair chosen depends on the boundary 
conditions. Outside the cylinders the asymptotic behavior determines which functions are 

(2) 

used. Since only the H m function has the behavior of an out-propagating cylindrical wave, 
the field of the scattered wave E sc has to be written using only#i 2) ofthe{H$>, H$} pair. 
Hence the scattered field is E sc = Xlm=-oo b m Hm (k o p)e m0 . Inside the cylinders we have to 
choose the { J m , Y m } pair because the Hankel functions are both singular at the origin. This 
behavior comes from the Y m part of the Hankel function (Hm = J m + iY m , Hm = Jm~iY m ) 
and thus only the J m functions can be kept for the interior. In this case the total internal 
field is E%* = Em=-oo a m Jm(k lP y me . 

After expressing the incident field in polar coordinates, i.e. in a Fourier-Bessel series, 
we can write the electric and magnetic field continuity conditions at the boundaries of each 
cylinder. This gives us a system of equations with the unknowns being the Fourier-Bessel 
coefficients of the scattered field outside the cylinders and total field inside the cylinders. The 
system can be simplified by using the relationship between the Bessel-Fourier coefficients 
of a field incident on a cylinder and the coefficients for the scattered and internal fields 

□ □□□ 



APPENDIX B: ERROR FUNCTION 



The optimization is based on the minimization of a functional, which is defined as the 
residual error between the calculated angular distribution of the normal component of the 
Poynting vector S and a distribution expressed as a target function. This error function 
is computed starting from the difference in intensity between the target and the result. 
In 2D the error is calculated along an observation line (often a circle around the group of 
cylinders). This line is divided into small portions and the normal component of the S vector 
is calculated in the center of each segment. 

Let the target function T(a) be the angular distribution of intensity exiting the circular 
observation region and S n (a) the normal component of the real part of the Poynting vector 
(S n (a) = S(a) • n(a) where n is the normal unit vector). In this space of functions defined 
on [0°, 360°] and having real values we can define a "distance" D between result S n (a) and 
target T(a) as D = ^ \S n (a) — T(a)[ l da. To properly evaluate the difference between 
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the target and result the functions T and S and have to be similarly normalized. 

In general, the exponent 7 can take any value. Choosing 7 = 1 assures that each im- 
provement is considered with the same weight. When choosing 7 > 1 improvements made 
in regions where the target and the results are very different influence the integral D more 
than a few smaller improvements in other regions. This means that 7 = 1 tends to ensure an 
uniform convergence while 7 > 1 favors reduction of major differences between target and 
result. The greater the numerical value of 7, the more important this effect becomes, while 
for < 7 < 1 the effect is reversed. And finally, a negative exponent 7 tends to push the 
solution further away from the target, in a manner which depends on the numerical value 
of 7. 

One could use different exponents for different stages of the iterative process. For example 
7 = 1 could be used at the beginning of the convergence procedure to avoid local minima. 
Later, the value of 7 could be increased to accelerate convergence towards a minimum. 
Furthermore, if one decides that this minimum is not sufficiently close to the target function, 
application of a negative exponent would divert the iterations from this local minimum to 
some intermediate point were 7 > 1 iterations could be restarted in the search for a better 
minimum. 

We use the real part of the Poynting vector in our calculation of S n (a). For most scat- 
tering directions this is the Poynting vector of the total field, with the exception of the 
input beam region where we use the scattered field in the first example (the top-hat target 
function) and the difference between the Poynting vector of the input field and total field in 
the second example (the cos 2 target function). 

The distances, or errors, D are computed numerically by dividing the observation circle 
into very small and equal portions and the integral replaced by a sum over these portions 



n* — v^ jv p 



Np Snjaj) _ Tjoj) 

N s N t 



where N s and N t are normalization factors. Different normal- 
ization methods can be used such as normalization to the maximum value, normalization 
to the sum of all values, or normalization to the sum of the squares. Since in our case the 
functions involve intensity we have chosen the normalization to be the sum i.e. the total 
power. 

Although not used in this paper, it is worth mentioning that when optimizing for a 
modal shape distribution (both amplitude and phase) a better choice for error function 
would be the overlap integral between the actual field distribution and the desired modal 
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field distribution. In this situation the target would be maximized. 
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FIG. 1: Starting configuration for the top-hat target function example. 56 dielectric rods with 
diameter d = 0.4 //m, index n = 1.5, in air. The incident beam, which propagates along the x axis 
(left to right), is a TM polarized gaussian beam with beamwidth 2a = 4 /xm, and wavelength A = 1 
jum. (a) Positions of the cylinders and distribution of the Poynting vector field. The observation 
circle, has a radius of 7 ^m and it is represented without the target window between 30 and 60 
degrees, (b) Computed angular intensity distribution (continuous line) at radius 7 ^m compared 
with the target distribution (dashed line). 



FIG. 2: Optimized configuration for the top-hat target function example. 56 dielectric rods with 
diameter, d = 0.4 ^m, index n = 1.5, in air. TM polarization gaussian beam with incident 
beamwidth 2a = 4 ^m, and wavelength A = 1 ^m. (a) Positions of the cylinders and distribution 
of the Poynting vector field after 9,700 iterations. The observation circle, has a radius of 7 /im and 
it is represented without the target window between 30 and 60 degrees, (b) Computed angular 
intensity distribution after 9,700 iterations (continuous line) at radius 7 /im compared with the 
target distribution (dashed line), (c) Contour plot of the electric field magnitude in relative units 
after 9,700 iterations, (d) Evolution of the error. 
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FIG. 3: Comparison with a periodic structure (PC) with 56 dielectric rods (n = 1.5) in air. 
The lattice constant and the angle are chosen so that for the incident wavelength A = l//m the 
Bragg diffraction condition is satisfied for 45°. (a) Positions of the cylinders and distribution of 
the Poynting vector field. The incident beam is propagating along the x axis (left to right). The 
observation circle, has a radius of 7 fim and it is represented without the target window between 
30 and 60 degrees, (b) Computed angular intensity distribution (continuous line) at radius 7 fim 
compared with the target distribution (dashed line). 



FIG. 4: Initial configuration for the cosine squared target function example. Si (n = 3.5) with 
26 cylindrical holes filled with Si02, TM gaussian beam incidence, beamwidth 2a = 1.5 fim, and 
wavelength A = 1.5 fim. (a) Positions of the cylinders and distribution of the Poynting vector 
field. The incident beam is propagating along the x axis (left to right). The observation circle, 
has a radius of 6 fim and it is represented without the target window between 30 and 60 degrees, 
(b) Computed angular intensity distribution (continuous line) at radius 6 ^m compared with the 
target distribution (dashed line). 



FIG. 5: Optimized configuration for the cosine squared target function example. Si (n = 3.5) 
having 26 cylindrical holes filled with Si02, TM gaussian beam incidence, beamwidth 2a = 1.5 
fim, and wavelength A = 1.5 fim. (a) Positions of the cylinders and distribution of the Poynting 
vector field. The incident beam is propagating along the x axis (left to right) after 7,000 iterations. 
The observation circle, has a radius of 6 fim and it is represented without the target window 
between 30 and 60 degrees, (b) Computed angular intensity distribution after 7,000 iterations 
(continuous line) at radius 6 fim compared with the target distribution (dashed line), (c) Contour 
plot of the magnetic field magnitude in relative units after 7,000 iterations, (d) Evolution of the 
error as a fucntion of iteration number. 



14 



FIG. 6: Comparison between 3D and 2D simulations, (a) The 3D structure simulated and the 
modal field of the input ridge waveguide. The thickness of the Si slab in the 3D simulation is 
0.6 fim and the (horizontal) mode size of the ridge waveguide is approximately w = 1.5 /im, the 
diameter of the holes is d = 0.4 jum. Incident light is TE polarized and has wavelength A = 1.5/im. 
(b) 3D Simulation results showing the Poynting vectors field in the middle horizontal section of the 
Si slab, (c) Comparison between the angular intensity distribution for the 3D and 2D simulations, 
(d) 2D simulation results, the Poynting vectors field. 



FIG. 7: Sensitivity analysis with respect to the frequency of the incident light for the optimized 
structure from Fig. |SJ Percentage change in error when the frequency of the light is modified from 
the original value /o = 200 THz . 



FIG. 8: Sensitivity analysis with respect to position for the optimized structure from Fig. [HJ The 
10 most sensitive cylinders are randomly moved along the x and y axis by a fixed length step 
(horizontal axis 10, 20, 30, 40 and 50 nm). The errors for 10 different randomly generated moved 
distributions are plotted for each step size. 
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